
#################################################################################################

#  Description: Replicating main analysis from Fleming, Thomas G. (Forthcoming) 
#               'Why Change a Winning Team? Explaining Post-Election Cabinet Reshuffles in Four
#                Westminster Democracies', Political Studies.
#  Author:      Thomas G. Fleming
#  Date:        10/09/2021
#  Note:        Requires file "reshuffles_data.csv" to be saved in same folder

#################################################################################################

rm(list=ls())

library(stargazer)  # for stargazer()
library(ggplot2)    # for ggplot()
library(MASS)       # for glm.nb()
library(ggpubr)     # for ggscatter()
library(DescTools)  # for PseudoR2()
library(margins)    # for margins()

library(rstudioapi) # sets location as wd
current_path <- getActiveDocumentContext()$path
setwd(dirname(current_path))

####################################
# Loads data and formats variables #
####################################

  data <- read.csv("./reshuffles_data.csv")
  
  data$election.date <- as.Date(data$election.date, format="%Y-%m-%d")  # formats as date
  data$cabinet.date <- as.Date(data$cabinet.date, format="%Y-%m-%d")    # formats as date
  
  data$selection.rule <- relevel(data$selection.rule, ref = "pm")       # sets reference category
  
################################################################################################
# Creates dependent variable - number of pre-election ministers moved within or out of cabinet #
################################################################################################

  data$cabinet.changes <- data$moved + data$left
    
####################################################################################
# Checking dispersion of variable for choosing count model - clearly overdispersed #
####################################################################################

  mean(data$cabinet.changes)
  var(data$cabinet.changes)
      
################################################################
# Table 1: Producing descriptive statistics for main variables #
################################################################
  
  varlist <- data[c("cabinet.changes", "seatshare", "seatshare.change", "preelection.total", "unavailable", 
                    "median.post.tenure", "median.cabinet.tenure", "new.PM", "selection.rule")]
  
  varlist$selection.rule <- ifelse(varlist$selection.rule == "caucus", 1, 0)

  lapply(varlist, mean, na.rm = TRUE)
  lapply(varlist, sd, na.rm = TRUE)
  lapply(varlist, min, na.rm = TRUE)
  lapply(varlist, max, na.rm = TRUE)
  
#####################################################################################
# Figure 1: Plots the bivariate relationship between cabinet changes and seat share #
#####################################################################################
  
  # seat share
  
  seatshare.fig <- ggscatter(data, x = "seatshare", y = "cabinet.changes", 
                             size = 0.5,
                             add = "reg.line", conf.int = TRUE, fullrange = TRUE)+
    theme_classic()+
    theme(axis.title.x = element_text(face = "bold"),
          axis.title.y = element_text(face = "bold"),
          panel.spacing.x = unit(4, "mm",),
          panel.border = element_rect(colour = "black", fill = NA, size = 0.7),
          axis.line = element_blank())+
    labs(y = "Cabinet changes\n",
         x = "\nRe-elected governing party seat share (%)")+
    facet_wrap(. ~ country, ncol = 4, scales = "free")
  
  # seat share change
  
  seatshare.change.fig <- ggscatter(data, x = "seatshare.change", y = "cabinet.changes", 
                                    size = 0.5,
                                    add = "reg.line", conf.int = TRUE, fullrange = TRUE)+
    theme_classic()+
    theme(axis.title.x = element_text(face = "bold"),
          axis.title.y = element_text(face = "bold"),
          panel.spacing.x = unit(4, "mm",),
          panel.border = element_rect(colour = "black", fill = NA, size = 0.7),
          axis.line = element_blank())+
    labs(y = "Cabinet changes\n",
         x = "\nChange in re-elected governing party seat share (%)")+
    facet_wrap(. ~ country, ncol = 4, scales = "free")
  
  # combining both seat share and seat share change into one figure
  
  bivariate.fig <- ggarrange(seatshare.fig, seatshare.change.fig, ncol = 1)
  ggsave("./fig1_bivariate.png", 
         width = 12, 
         height = 5.5, 
         units = "in")
  
###################################################################################################
# Table 2: Models the number of new appointments as a function of seat share and various controls #
###################################################################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.bivariate <- glm.nb(cabinet.changes ~
                                  seatshare +
                                  preelection.total,
                                data = data)
  seatshare.bivariate.R2 <- PseudoR2(seatshare.bivariate, which = "McFadden")
  
  # adding various control variables
  
  seatshare.controls <- glm.nb(cabinet.changes ~
                                 seatshare +
                                 preelection.total +
                                 unavailable +
                                 median.post.tenure +
                                 median.cabinet.tenure +
                                 new.PM +
                                 selection.rule,
                               data = data)
  seatshare.controls.R2 <- PseudoR2(seatshare.controls, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.countryFE <- glm.nb(cabinet.changes ~
                                  seatshare +
                                  preelection.total +
                                  unavailable +
                                  median.post.tenure +
                                  median.cabinet.tenure +
                                  new.PM +
                                  selection.rule +
                                  as.factor(country),
                                data = data)
  seatshare.countryFE.R2 <- PseudoR2(seatshare.countryFE, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.partyFE <- glm.nb(cabinet.changes ~
                                seatshare +
                                preelection.total +
                                unavailable +
                                median.post.tenure +
                                median.cabinet.tenure +
                                new.PM +
                                selection.rule +
                                as.factor(party),
                              data = data)
  seatshare.partyFE.R2 <- PseudoR2(seatshare.partyFE, which = "McFadden")
  
  # creating combined table of coefficients
  
  vars = c("seatshare", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", 
           "selection.rulecaucus", "Constant")
  
  stargazer(seatshare.bivariate, 
            seatshare.controls, 
            seatshare.countryFE, 
            seatshare.partyFE,
            type = "html",
            title = "Table 2: Negative binomial regression of post-election cabinet changes",
            out = "./table2.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Seat share", "Total ministers", "Unavailable ministers", "Median tenure (current post)", "Median tenure (any post)",
                                 "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(1)", "(2)", "(3)", "(4)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.bivariate), nobs(seatshare.controls), 
                               nobs(seatshare.countryFE), nobs(seatshare.partyFE)),
                             c("Pseudo R2", round(seatshare.bivariate.R2, digits = 3),
                               round(seatshare.controls.R2, digits = 3),
                               round(seatshare.countryFE.R2, digits = 3),
                               round(seatshare.partyFE.R2, digits = 3))),
            font.size = "small")
  
################################################################################################
# Table 3: Repeats previous models with *change in* seat share as the new independent variable #
################################################################################################
  
  # bivariate (controlling just for total ministers)
  
  seatshare.change.bivariate <- glm.nb(cabinet.changes ~
                                         seatshare.change +
                                         preelection.total,
                                       data = data)
  seatshare.change.bivariate.R2 <- PseudoR2(seatshare.change.bivariate, which = "McFadden")
  
  # adding various control variables
  
  seatshare.change.controls <- glm.nb(cabinet.changes ~
                                        seatshare.change +
                                        preelection.total +
                                        unavailable +
                                        median.post.tenure +
                                        median.cabinet.tenure +
                                        new.PM +
                                        selection.rule,
                                      data = data)
  seatshare.change.controls.R2 <- PseudoR2(seatshare.change.controls, which = "McFadden")
  
  # adding country fixed effects
  
  seatshare.change.countryFE <- glm.nb(cabinet.changes ~
                                         seatshare.change +
                                         preelection.total +
                                         unavailable +
                                         median.post.tenure +
                                         median.cabinet.tenure +
                                         new.PM +
                                         selection.rule +
                                         as.factor(country),
                                       data = data)
  seatshare.change.countryFE.R2 <- PseudoR2(seatshare.change.countryFE, which = "McFadden")
  
  # replacing country fixed effects with party fixed effects
  
  seatshare.change.partyFE <- glm.nb(cabinet.changes ~
                                       seatshare.change +
                                       preelection.total +
                                       unavailable +
                                       median.post.tenure +
                                       median.cabinet.tenure +
                                       new.PM +
                                       selection.rule +
                                       as.factor(party),
                                     data = data)
  seatshare.change.partyFE.R2 <- PseudoR2(seatshare.change.partyFE, which = "McFadden")
  
  # combining results in single regression table
  
  vars = c("seatshare.change", "preelection.total", "unavailable", "median.post.tenure", "median.cabinet.tenure", "new.PM", 
           "selection.rulecaucus", "Constant")
  
  stargazer(seatshare.change.bivariate, 
            seatshare.change.controls, 
            seatshare.change.countryFE, 
            seatshare.change.partyFE,
            type = "html",
            title = "Table 3: Negative binomial regression of post-election cabinet changes",
            out = "./table3.html",
            keep = paste0("^", vars, "$"),
            order = paste0("^", vars, "$"),
            covariate.labels = c("Seat share change", "Total ministers", "Unavailable ministers", "Median tenure (current post)", "Median tenure (any post)",
                                 "PM's first election", "Caucus selection", "Constant"),
            omit.stat = c("n", "ll", "theta", "aic"),
            no.space = TRUE,
            dep.var.caption = "DV: Ministers moved from pre-election post",
            dep.var.labels.include = FALSE,
            model.numbers = FALSE,
            column.labels = c("(5)", "(6)", "(7)", "(8)"),
            add.lines = list(c("Country FE", "No", "No", "Yes", "No"),
                             c("Party FE", "No", "No", "No", "Yes"),
                             c("Observations", nobs(seatshare.change.bivariate), nobs(seatshare.change.controls), 
                               nobs(seatshare.change.countryFE), nobs(seatshare.change.partyFE)),
                             c("Pseudo R2", round(seatshare.change.bivariate.R2, digits = 3),
                               round(seatshare.change.controls.R2, digits = 3),
                               round(seatshare.change.countryFE.R2, digits = 3),
                               round(seatshare.change.partyFE.R2, digits = 3))),
            font.size = "small")

######################################################
# Figure 2: Plots the marginal effects of seat share #
######################################################
    
  # Calculates marginal effects from seatshare models
  
  model1.ME <- summary(margins(seatshare.bivariate, variables = "seatshare"))
  model2.ME <- summary(margins(seatshare.controls, variables = "seatshare"))
  model3.ME <- summary(margins(seatshare.countryFE, variables = "seatshare"))
  model4.ME <- summary(margins(seatshare.partyFE, variables = "seatshare"))
  
  seatshare.MEs <- rbind(model1.ME, model2.ME)
  seatshare.MEs <- rbind(seatshare.MEs, model3.ME)
  seatshare.MEs <- rbind(seatshare.MEs, model4.ME)
  
  # Adds model names and sets order of levels for figure presentation
  
  seatshare.MEs$model <- c("Model 1", "Model 2\n(Controls)", "Model 3\n(Country FE)", "Model 4\n(Party FE)")
  
  seatshare.MEs$model <- factor(seatshare.MEs$model, 
                                levels = seatshare.MEs$model[order(seatshare.MEs$model)])
  
  # Adds 90% confidence intervals
  
  seatshare.MEs$upper90CI <- seatshare.MEs$AME + 1.645*seatshare.MEs$SE
  seatshare.MEs$lower90CI <- seatshare.MEs$AME - 1.645*seatshare.MEs$SE
  
  # Produces plot of these marginal effects
  
  ggplot(data = seatshare.MEs, aes(y = model, x = AME, xmin = lower90CI, xmax = upper90CI))+
    geom_pointrange()+
    scale_y_discrete(limits = rev(levels(seatshare.MEs$model)))+
    xlim(-0.45,0.2)+
    geom_vline(xintercept = 0, linetype = "dashed")+
    labs(y = NULL, x = "\nAverage Marginal Effect")+
    theme_classic()+
    theme(axis.title.x = element_text(face = "bold"),
          panel.border = element_rect(colour = "black", fill = NA, size = 0.7), # adds rectangle
          axis.line = element_blank())                                          # removes axis lines because they reinforce the rectangle 
  
  ggsave("./fig2_seatshare_MEs.png", 
         width = 6, 
         height = 3,
         units = "in")
  
#############################################################
# Figure 3: Plots the marginal effects of seat share change #
#############################################################
  
  # Calculates marginal effects from seatshare change models
  
  model5.ME <- summary(margins(seatshare.change.bivariate, variables = "seatshare.change"))
  model6.ME <- summary(margins(seatshare.change.controls, variables = "seatshare.change"))
  model7.ME <- summary(margins(seatshare.change.countryFE, variables = "seatshare.change"))
  model8.ME <- summary(margins(seatshare.change.partyFE, variables = "seatshare.change"))
  
  seatshare.change.MEs <- rbind(model5.ME, model6.ME)
  seatshare.change.MEs <- rbind(seatshare.change.MEs, model7.ME)
  seatshare.change.MEs <- rbind(seatshare.change.MEs, model8.ME)
  
  # Adds model names and sets order of levels for figure presentation
  
  seatshare.change.MEs$model <- c("Model 5", "Model 6\n(Controls)", "Model 7\n(Country FE)", "Model 8\n(Party FE)")
  
  seatshare.change.MEs$model <- factor(seatshare.change.MEs$model, 
                                       levels = seatshare.change.MEs$model[order(seatshare.change.MEs$model)])
  
  # Adds 90% confidence intervals
  
  seatshare.change.MEs$upper90CI <- seatshare.change.MEs$AME + 1.645*seatshare.change.MEs$SE
  seatshare.change.MEs$lower90CI <- seatshare.change.MEs$AME - 1.645*seatshare.change.MEs$SE
  
  # Produces plot of these marginal effects
  
  ggplot(data = seatshare.change.MEs, aes(y = model, x = AME, xmin = lower90CI, xmax = upper90CI))+
    geom_pointrange()+
    scale_y_discrete(limits = rev(levels(seatshare.change.MEs$model)))+
    xlim(-0.2,0.2)+
    geom_vline(xintercept = 0, linetype = "dashed")+
    labs(y = NULL, x = "\nAverage Marginal Effect")+
    theme_classic()+
    theme(axis.title.x = element_text(face = "bold"),
          panel.border = element_rect(colour = "black", fill = NA, size = 0.7), # adds rectangle
          axis.line = element_blank())                                          # removes axis lines because they reinforce the rectangle 
  ggsave("./fig3_seatsharechange_MEs.png", 
         width = 6, 
         height = 3,
         units = "in")
  
########################################################################################## 
  
  rm(list=ls())
  